delta_u = zeros(N,2);
for i = 1:N-1
    delta_u(i,:) = u_piao(i + 1,:) - u_piao(i, :);
end

plot(delta_u(:,1));
% plot(delta_u(:,2));


